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FINITE ELEMENT ANALYSIS AND COMPUTER GRAPHICS VISUALIZATION 
OF FLOW AROUND PITCHING AND PLUNGING AIRFOILS 

By Theodore Bratanow* and Akin Ecer** 
University of Wisconsin-Milwaukee 


SUMMARY 


A general computational method for analyzing unsteady flow around 
pitching and plunging airfoils was developed. The finite element method 
was applied in developing an efficient numerical procedure for the solu- 
tion of equations describing the flow around airfoils. Numerical results 
were employed in conjunction with computer graphics techniques to produce 
visualization of the flow. The investigation involved mathematical model 
studies of flow in two phases: analysis of a potential flow formulation 

and analysis of an incompressible, unsteady, viscous flow from Navier- 
Stokes equations. 

For the analysis of the potential flow formulation of the problem, 
a finite element representation of Laplace's equation was applied. The 
potential flow model was extensively employed in reproducing patterns of 
unsteady flow characteristics, in the form of computer-drawn illustrations 
and computer-drawn movies. Using the finite element method again, mathe- 
matical models were further developed to analyze unsteady viscous flow 
around pitching and plunging airfoils from solutions of the vorticity trans- 
port equation. A system of first order nonlinear matrix differential equa- 
tions, resulting from the finite element formulation, was integrated with 
respect to time to obtain vorticity distribution around airfoils. A system 
of Lagrange multipliers was applied in representing boundary conditions 
around airfoils, which change with respect to time and space. Since time- 
dependent boundary conditions on airfoils were treated separately in the 
formulation, the applied method produced efficient solutions of the vorti- 
city transport equation, thus significantly reducing the computational 
effort . 

The results include computer -drawn movies and graphical representa- 
tions of streaklines, streamlines, and vorticity distributions around 
pitching and plunging airfoils. Various features of the investigation, 
including a detailed description of both mathematical models as well 
as applied computational and computer graphics techniques, are presented. 
Accuracy and advantages of the developed procedure are evaluated. 


* Professor, ** Assistant Professor, Department of Mechanics, 

College of Engineering and Applied Science. 


INTRODUCTION 


The general motion of a flexible wing may consist of pitching and 
plunging components . The flow around the wing can be of an unsteady and 
three-dimensional character. A complete and clear description of such 
three-dimensional flows involves analytic methods applied in unsteady 
aerodynamics. Solutions of three-dimensional mathematical representa- 
tions of such problems, as well as usable experimental data to support 
analytic work, are not easily obtained. Up-to-date analytic efforts 
have not been successful in producing desirable results. Difficulties 
have also been encountered in experimental work with scaling and measure- 
ing the required data. The analysis of unsteady aerodynamics problems 
involves solutions of nonlinear three-dimensional partial differential 
equations. Closed form solutions of such equations are difficult and 
exist only for a limited number of cases. Even with existing advanced 
computational facilities, numerical solutions of such problems require 
excessive computer time. 

The purpose of the investigation was to develop mathematical 
models for analysis and visualization of general unsteady flow problems 
around pitching and plunging airfoils. The numerical results, defining 
the motion of the airfoil and the flow, were used in conjunction with 
computer graphics techniques to produce adequate visualization in 
time and space of the flow around the oscillating airfoil. 

The overall mathematical procedure involves two major steps, which 
are closely related: 

I. Mathematical model analysis 

II. Generation of computer-drawn illustrations of flow patterns and 
computer-drawn movies for visualization of the flow around air- 
foils. 

The computational procedure for the potential and unsteady viscous flow 
analyses can be summarized as follows: 

defining the airfoil motion and flow parameters at each time step 
as input 

obtaining the velocity distribution around the airfoil at each 
time step from the mathematical models 

numerical integration of the obtained velocity distribution to 
describe desired characteristics of the flow 

illustration of the flow characteristics using computer graphics 
techniques . 

The analytical formulations were established in such a way that they 
were suitable for expansion in desired directions. Two-dimensional solu- 
tions of potential flow formulations proved to be most fruitful as a basis 
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for gaining experience and then gradually increasing the complexity of the 
mathematical model. The next step of the investigation was the analysis of 
incompressible, unsteady and viscous flow from Navier-Stokes equations. The 
analytically obtained flow patterns were then compared with patterns obtained 
by others and with actual physical flow. Nonlinearities were treated both 
analytically and computationally. A direct method was applied for efficient 
numerical solution of the equations representing the flow around the airfoil 
from the standpoint of computer time consumption. The finite element method 
was applied in a broad sense, beyond the limitations of the classical var- 
iational approach to solve the Navier-Stokes equations. The overall results 
were satisfactory. 

The presented results are from the analysis of two-dimensional flow 
around a pitching and plunging airfoil. The analytical tools are in a 
form capable of treating general flow problems. 


A 

A. 
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SYMBOLS 

2 2 

area of defined region, in. ( m ) 

2 2 

area of a triangular element, in. ( m ) 
symmetric matrix 

symmetric matrix for one triangular element 

symmetric matrix containing matrix A and the 
boundary conditions 

width of a typical triangular element in figure 4, 
in. ( m ) 


b 


rectangular matrix representing the boundary con- 
ditions on the airfoil 


^ row vector for the entire grid 

row vector for a triangular element 

c_ rectangular matrix specifying the position of the 

nodes on the free stream boundary 

^ column vector specifying the stream functions for 

the free stream boundary 
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0 
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g( x>y ) ■ a function specifying the boundary con- 
ditions 

rectangular matrix containing matrix = C ^ 2. ^ 

height of a typical triangular element in figure 4, 
in. ( m ) 

column vector containing vector ~ ( 0 2 ^ 

identity matrix 
null matrix 

2 

local static pressure, psi ( N/m ) 


£ 

Q 
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column vector, £ = A 


-1 


*-l 

rectangular matrix, £ = A g^ 

A 2 2-2 -2 2 
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Q( x,y ), Q = - pi — 2~ * — 2~ 

\3x 9x9y 9y 


column vector representing a discrete finite element 
formulation of Q 

a column vector 
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a column vector 
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x,y,z 


X. . , y . . 

ij ij 


X,Y 


u,v 


u 
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time, sec. 

coordinates, in. ( m ) 
coordinate differences, in. C ) 

body forces in x and y directions respectively, 
Ibf ( N ) 

velocity components of the fluid in the x and y 
directions respectively, in. /sec. ( m/sec. ) 

velocity vector, in. /sec. ( m/sec. ) 

3 3 

del operator, V = — + 
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6 
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p 

3 > 

i 


( 0 . 

—1 


II * o2 3 3 

Laplace operator, V = — ^ ^ 

3x^ 

an eigenvalue of the triangular finite element grid 
column matrix representing the Lagrangian multipliers 
largest eigenvalue of a triangular element 
eigenvalue of a triangular element 

2 2 

kinematic viscosity coefficient, in. /sec. (m /sec.) 

3 3 

density of air, Ibra/in. (Kg/m ) 

4i(ip) - quadratic function 

2 2 

>^(x,y) - stream function, in. /sec. (m /sec.) 

column vector representing stream function values 

2 2 

at grid intersection points, in. /sec. (M /sec.) 
column vector, 4^* = (ij^ ^) 

2 

values of the stream function at a boundary, in. /sec. 
(m^/sec. ) 

column vector representing stream function values at 

2 

the corners of a triangular element, in. /sec. 

(m^/sec . ) 

<jj(x,y) - vorticity function (sec. ^) 

column vector representing vorticity values at the 
nodal points (sec. ^) 

column vector representing vorticity values at the 
corners of a triangular element (sec. ^) 
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Subscripts 

b 

c 

h 

i 

3 

max 

o 


P 

q 

s 

t 

V 

x,y 

Superscripts 

n 

t 


at the airfoil boundary 

at the free stream boundary 

horizontal 

i-th triangle 

j-th triangle 

maximum 

at time t 

o 

p-th node 
q-th node 

at the airfoil and free stream boundaries 

at time t 

vertical 

derivative with respect to x or y 

at the n-th time step 
transpose 
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COMPUTER ANALYSIS AND VISUALIZATION OF FLOW PROBLEMS 


Background 

A complete three-dimensional analysis of the flow around a pitching 
and plunging wing in flight would demonstrate the complexity of problems 
involving unsteady aerodynamics. Applications of numerical methods in 
solving unsteady aerodynamics flow problems have received considerable 
attention during recent years [refs. 1, 2, 3, 4, 5, 6]. Attempts have 
also been reported of solving three-dimensional flow problems numerically 
by considering them as two-dimensional ones and then including only lin- 
ear effects of the third dimension [refs. 7, 8, 9]. Some approximate 
solutions have been obtained from linearized forms of three-dimensional 
flow formulations. In each instance, however, a clear understanding of 
the kinematics of the three-dimensional flow was required in order to 
obtain realistic results. Some investigations have involved application 
of finite difference techniques and have yielded encouraging but limited 
results. Applications of finite difference techniques did not become very 
popular because of the amount of computer time required for satisfactory 
results. Also, in some reported investigations theoretical and exper- 
imental methods have been used conjunctional ly. 

The finite element method was originally applied in the field of 
solid mechanics. Availability of variational formulations for some of 
the basic solid mechanics problems encouraged researchers to take advan- 
tage of such formulations in applying the finite element method. It 
should be pointed out, however, that variational formulations may not be 
the only basis for finite element applications. The application is essen- 
tially a representation of a continuum by discrete finite elements and 
setting up governing equations of the dynamics problem in terms of discrete 
displacements, velocities and accelerations at elemental nodal points. 

The finite element method has already been applied for solution of 
flow problems [refs. 10, 11, 12, 13, 14]. Early applications favored 
problems for which variational formulations existed. Examples are: 
seepage flow, heat conduction and potential flow problems which involve 
solutions of Laplace's equation. Many reported attempts at solving gen- 
eral flow problems encountered difficulties with the nonlinearities in- 
volved in the Navier-Stokes equations. Some flow problem formulations 
[refs. 15, 13] have been based on the concept of "local potential", as 
introduced by Glansdorff and Prigogine [ref. 16] . In such applications 
the energy and momentum can be minimized in a finite region, for which 
equations of motion can be obtained. 

Oden [ref. 12] has applied the finite element method for general prob- 
lems involving slow viscous flow. The variational formulation was based on 
minimum viscous dissipation as described by the Helmholtz equation. The 
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solutions are of the type obtained from Dirichlet's problem. The formula- 
tion is concise, and accuracy of the results is satisfactory. Bowley and 
Prince [ref. 10] developed a numerical procedure which uses a finite ele- 
ment technique. They applied the integral form of the time dependent equa- 
tions of mass, (x- and y-) momenta, and energy conservation to various 
subdomains of the whole region. The results were also satisfactory. 

In 1963 Fromm [ref. 3] applied the Helmholtz vorticity transport equa- 
tion in conjunction with the finite difference method to analyze unsteady 
incompressible flow at various Reynolds numbers around stationary objects. 
Derivatives in the Helmholtz equation were replaced with central differ- 
ences both in space and time. A successive approximation technique was 
applied to obtain the streamfunction, which satisfies Poisson's equation. 
Fromm's results demonstrated in a spectacular way the feasibility of such 
a technique; however, the applied numerical solution contained already 
known deficiencies of the finite difference method. The applied iterative 
scheme slows down the computational process. Even more significant is the 
difficulty connected with the representation of a moving obstacle having a 
curved boundary, as is the case with an oscillating airfoil. Obtaining a 
reasonably accurate solution of such problems by applying finite difference 
methods seems to be uneconomical. 

Various experimental techniques have been applied in attempts to vis- 
ualize and measure unsteady flow around airfoils in the laboratory; only 
limited results were obtained and correlation was poor. It is not known 
if other efforts have been made for the flow visualization problem of a 
pitching and plunging airfoil in conjunction with a theoretical analysis, 
involving a finite element method and computer graphics techniques. 


Mathematical Analysis 

The procedure followed in the mathematical model analysis involves 
mainly: 


solutions of a potential flow representation of the flow 
patterns around oscillating airfoils 

solutions of Navier-Stokes equations for representation of 
patterns of unsteady, incompressible, viscous flow around 
pitching and plunging airfoils. 

Potential flow analysis . - Numerical techniques were applied to ob- 
tain solutions of the potential flow representation. The procedure for 
the potential flow analysis will now be summarized. Streamfunctions in 
the flow field were defined as 


8 



dip 

dy 


V 


( 1 ) 


Sx 

where the streamfunctions satisfy Laplace's equation 

V^<p = 0 (2) 

for prescribed conditions 


^5 = g (3) 

on airfoil boundaries and on imaginary boundary lines representing the 
undisturbed flow. 

The developed numerical procedure for solution of the system of equa- 
tions (1), (2), and (3) was based on the finite element method [ref. 17]. 
Details of the mathematical model and the numerical analysis are summari- 
zed in Appendix A. 

Unsteady, viscous flow analysis . - Solutions for unsteady, incompress- 
ible and viscous flow are obtained from Navier-Stokes equations as a re- 
finement of the potential flow model. The system of vorticity transport 
equation and Poisson's equation applied in the analysis of flow around an 
oscillating airfoil is 


3u) ^ 3 

3t 3x 





0 


(4) 


V^ip = - 0) (5) 

Equations (4) and (5) were solved simultaneously for the same time step, 
whereby equation (4) was integrated with respect to time, to obtain 
vorticity distribution around oscillating airfoils. The finite element 
method was applied again in obtaining a numerical solution. The devel- 
oped numerical method is quite general and can be applied in analyzing 
flow around obstacles of any shape under arbitrary motion. Detailed 
description of the mathematical model as well as the applied numerical 
method is presented in Appendix B. 


Details of the Mathematical Analysis 

Figures la and lb illustrate the important features of the applied 
analysis of mathematical models for potential flow and unsteady viscous 
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flow. The two- and three-dimensional potential flow formulations and 
analyses are very useful first approximations of flow patterns around 
oscillating airfoils [ref. 11]. However, the representations of the 
flow are general and do not account for the effects of the boundary 
layers . 

The finite element method was applied for both the numerical solution 
of the potential flow formulation, equations (1), (2), (3), and the 
numerical solution of the unsteady viscous flow formulation, equations (4) 
and (5). The numerical procedure for the treatment of Poisson's equation, 
equation (5) , for both potential and unsteady viscous flow is direct and 
requires no iteration. The input data for the mathematical models 
describe the geometry of the obstacle and the free stream characteristics 
of the flow in a most general form. The exact shape of an airfoil and 
its motion are represented in their actual form, without any approximation. 

Initially, the vorticity distribution for each of the cases was 
unknown. Therefore, at the beginning of each numerical integration pro- 
cess, the vorticity distribution was assumed to be unity at a point 
approximately one half chord length upstream of the airfoil and zero for the 
remaining region around the airfoil. With the initially assumed vorticity 
distribution, Poisson's equation was solved for stream, functions and then 
the vorticity transport equation was integrated with respect to time. 

The numerical integration was advanced with respect to time and the flow 
characteristics at each time step were recorded. This process was con- 
tinued until the numerical integration reached a steady state; i.e., 
when the changes in the vorticity distribution with respect to time were no 
longer significant. The steady state condition was then used as an initial 
condition ' for analyzing the unsteady flow around oscillating airfoils. 

The vorticity distribution around the airfoil was obtained from the 
numerical integration with respect to time of a system of first order 
nonlinear matrix differential equations. The time-dependent boundary 
conditions in the solution of the vorticity transport equation and 
Poisson's equation were treated in a way to minimize computational cost. 

The boundary conditions around the obstacle in the flow field, which vary 
in time and space coordinates, were represented by resorting to the 
Lagrange multiplier technique. Doing so allowed for a separation of the 
boundary conditions in both' the solution of the vorticity transport equation 
and Poisson's equation. Terms which are independent of time were kept in 
computer storage throughout the entire numerical integration process in 
order to avoid excessive computer time consiomption. 

For the solution of the potential flow, the region on and around the 
obstacle in the flow field was defined by a finite element triangular grid 
as shown in figures 2a, 2b, and 2c. For each of the triangular elements 
in this grid, the variational form of Laplace's equation was developed. 
Around the leading edge of the airfoil, where the variation of the stream 
functions is most intense, a finer mesh was adopted. As the airfoil under- 
goes pitching and plunging motion, the entire grid can be rotated by simply 
changing the boundary conditions on the outside rectangle representing the 
free stream flow. 
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Another mesh, shown in figure 3, was used for the unsteady viscous 
flow analysis obtained from the Navier-Stokes equations. As can be seen 
from the figure, the arrangement of this mesh is not dependent on the 
positions of the airfoil. The boundary conditions on the airfoil are 
related to this mesh by simply setting prescribed values of stream func- 
tions and vorticities at certain points of the grid, using Lagrange 
multipliers. In the analysis of the flow around a circle, for example, 
the geometry of the circular obstacle is superimposed on the mesh. The 
motion of the obstacle can be described by redefining the boundary 
conditions on the obstacle for equations (4) and (5), at each time step of 
the numerical integration of equation (4) . These boundary conditions are 
then applied on the system of matrix differential equations, as obtained 
from the finite element representation of the governing differential 
equation. Solutions of the discretized form of Poisson's equation and 
the vorticity transport equation, at each time step, are then reduced to 
the solution of a system of linear algebraic equations. 

Lising a Lagrange multiplier technique the system of linear algebraic 
equations can be written in matrix form as 


A 

1 


s 


S 

b 

0 


6 


£ 


As the numerical integration proceeds, matrices ^ and ^ change with 
respect to time, depending on the instantaneous position of the obstacle. 
However, matrix A does not change with respect to time and can be inverted 
right at the beginning of the analysis and stored in the computer. The 
solution vector £ represents the streamfunctions in the solution of 
Poisson's equation and the rate of change of vorticities in the solution 
of the vorticity transport equation. 

The solution then consists of; 

assembling matrix which represents the vorticity distribution 
for the solution of Poisson's equation and the unsteady and vis- 
cous terms of the vorticity transport equation 

defining matrix b for a given shape of the obstacle and its 
instantaneous position 

. performing numerical integration of the discretized form of the 

vorticity transport equation. 

In the case of steady flow, the problem is simplified, since matrix 
^ does not change with respect to time. Considering that the order of ma- 
trix A turns out to be much larger than that of matrix the presented 

technique yields satisfactory results. 
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The Computer-Generated Movies 

Computer-generated movies offer definite advantages in visualiza- 
tion of solutions of time-dependent engineering problems. An exact 
interpretation of time-dependent results is faciliated by such visual- 
ization as these movies provide. Streamlines, streaklines, and pathlines 
were used in the flow visualization studies. Although such kinematic 
concepts are effective in flow visualization, a more forceful representa- 
tion of flow patterns is a movie showing the actual flow as it occurs; 
details of the motion can thus be depicted which are not easy to observe 
directly. 

Computer- generated movies are a direct product of analysis of 
mathematical models, for which assumptions and problem conditions are well 
defined. The accuracy of results obtained using computer graphics 
visualization techniques depends only on the accuracy of the developed 
mathematical model. The various parameters defining the motion of the 
fluid can be visualized in time and space. The observation can take place 
in any reduced time system, and details of the flow can be magnified and 
projected in any desired coordinate system. Simply by adjusting input 
parameters in the computer program, one can obtain a desirable view of 
the process and observe the effects of the individual physical parameters 
of the mathematical model. Corrections and improvements in the mathe- 
matical model can be observed immediately from the computer -drawn movies. 

Once the analysis of the mathematical model for the motion of the 
oscillating airfoil and the flow renders acceptable results, these results 
are used as input to produce computer-drawn movies. The computer graphics 
procedures developed for the flow visualization analysis were summarized 
in two main steps, as shown in figure Ic. The triangular grids shown in 
figures 2a, 2b, and 2c, which were used for the finite element analysis of 
potential flow, were used again for determining the motion of the fluid 
particles. The computer program first records the position of a fluid 
particle, which starts from a reference position, by determining in which 
triangular element the particle is located. The magnitude and direction 
of the velocity at that point is obtained from the output of the mathematical 
model; the motion of the particle is defined through numerical integration 
of this velocity vector with respect to time. Then, using analytic and 
computer graphics techniques the positions of a large number of fluid 
particles and the position of the obstacle at a particular respective time 
instant are projected on a frame. 

Results can be compared and correlated with actual flow information 
obtained from experimental studies and the mathematical model can again 
be improved. Such improvements are not easily made in an experimental 
procedure. Thus, a direct link can be established between experimental 
and theoretical studies. 

To achieve a satisfactory yet economical representation of the problem, 
it was necessary to scrutinize each individual sequence of the two analyses 
and their interrelation. A 4020 Stromberg-Carlson data graphics device was 
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used to generate the computer -drawn movies. The major portion of the 
computational work was carried out at the University of Wisconsin-Milwaukee 
on a UN I VAC 1108 computer. 

Discussion of the Results 

The results presented in this report were obtained from two-dimensional 
formulations of flow around airfoils in conjunction with computer graphics 
techniques. Physical and numerical aspects of unsteady flow problems 
were investigated. The presented data were chosen to show the suitability 
of the developed method for visualization of flow around pitching and 
plunging airfoils. The actual form of the airfoil was accurately represent- 
ed. The airfoil used in the analysis for flow visualization was NACA 0012 
airfoil. Several cases involving application of Poisson's equation and 
Navier-Stokes equations were analyzed. The results include streamlines, 
streaklines, vorticity distributions, and computer-generated movies. 

The movies depict accurately the effects of the pitching and plunging 
motion of the airfoil on the overall flow patterns around the airfoil. The 
production cost for the sample movies was quite moderate. The longest 
movie consisted of six sections, each 15 seconds long, and depicted the 
following cases: 

1) pitching of the airfoil through angles of attack ±8° 

2) plunging of the airfoil through ±85% of the airfoil chord thickness 

3) pitching of ±8° and plunging of ±85% chord thickness (the 

basic case) 

3a) same as in case 3 except that the plunging is ±42% of chord 
thickness 

3b) same as in case 3 except that the free stream velocity is 
reduced by 50% 

3c) same as in case 3 except that the pitching and plunging motion 
is 180° out of phase. 

An illustration of a few sample frames from a movie are shown in figure 4. 

The vorticity distribution around a NACA 0012 airfoil was obtained from 

the solution of the vorticity transport equation and Poisson's equation, as 
a direct output of the computer program. The different physical conditions 
for the numerical integration of the equations were; 

Reynolds numbers 10^, 10^, 10^ 

pitching amplitude 0° (static), 0°±4°, 0°±8° 

plunging amplitude ±42% and ±85% of chord thickness 
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different initial location of the starting vortex 

. timewise variation of the flow parameters. 

Computer graphics techniques were used to produce visualization of 
the vorticity distribution around the airfoil at successive time steps. 

The numerical results for the vorticity distribution around the airfoil, 
at various time steps, were plotted using a contour routine, as shown in 
figures 5-11. As can be seen from these figures, the initially imposed 
vortex diffuses around the airfoil with advancing time. The boundary 
conditions around a moving airfoil were described quite accurately during 
the calculations as can be seen for example from figure 11. The moving 
vortex in this figure follows quite closely the motion of the airfoil. 

The sample cases in figures 5-11 indicate: 

. difference in the overall behavior of the flow at different 
Reynolds numbers and angles of attack 

. influence on flow patterns due to pitching and plunging os- 
cillations of the airfoil 

I 

. behavior of a vortex after it strikes an airfoil 

. growth and decay of unsteadiness of the flow in areas where 
vorticity distribution increases. 

Figures 12 and 13 show computer-drawn illustrations of streamlines and 
streaklines. The streamlines were obtained from calculated values of the 
stream functions at particular time intervals by simply drawing contour 
lines. In figure 12 unsteady flow streamlines obtained from solutions of 
Navier-Stokes equations are shown at various angles of attack of the airfoil. 
Streaklines, as shown in figure 13, were obtained by plotting instantaneous 
positions of a numer of flow particles, which have passed through a fixed 
point in the free stream. After the velocity distribution around the air- 
foil at each time step was determined these values were integrated to obtain 
the position of each fluid particle. Comparisons of computer time required 
for the application developed here and finite difference applications have 
demonstrated the advantage of the developed procedure (See table 1). 


Concluding Remarks 

The developed computational method has dealt effectively with dif- 
ficulties involved in solutions of equations representing unsteady flow 
patterns around pitching and plunging airfoils. The cqmbined application 
of finite element method and the chosen numerical procedure, for both the 
analytic and graphics calculations, have provided a strong and flexible 
tool for producing desirable results at a reasonable computational cost. 
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The important advantages of the developed computational method can 
be summarized as follows: 

solution of Poisson's equation by a direct method, without any 
iteration, for both potential and unsteady viscous flow analyses 

. separate treatment of time -dependent boundary conditions related 
to oscillating airfoils, for both solution of Poisson's equation 
and numerical integration of the vorticity transport equation 

. capability to relate mathematical flow analysis directly to flow 
visualization 

. capability to represent boundary conditions on an obstacle of 
arbitrary shape. 

Experimental results applicable to flow visualization studies in 
unsteady aerodynamics are limited. The developed procedure for flow 
visualization has the potential of linking more effectively theoretical and 
experimental studies by filling the wide gap, which exists in many 
instances between results of the two investigative methods. The unsteady 
flow problem around an oscillating rotor blade airfoil is a most 
instructive example for such applications. A complete package of analysis 
and visualization of the flow around a pitching and plunging airfoil offers 
a new tool in aerodynamic studies. 
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APPENDIX A 


FINITE ELEMENT ANALYSIS OF POTENTIAL FLOW 
AROUND AN OSCILLATING AIRFOIL 


Method of Solution 

The solution of the problem represented by equations (2) and (3) can 
be reduced to finding the minimum of the integral 



(Al) 


where is a function of x and y which satisfies the boundary conditions 
in equation (3) [refs. 18, 19, 20]. 


The finite element method was applied for evaluation of the integral 
in equation (Al). The volume of air around the airfoil and the airfoil, 
accordingly, were represented by a series of triangles as shown in figure 
2 . 


The curved geometry of the airfoil boundary was considered in assem- 
oling the finite element grid. Once the finite element grid was devel- 
oped, it was assumed that the values of the streamfunctions vary line- 
arly over the area of each of the finite element triangles. Although a 
finite element model of a higher order could be used, in this case even a 
lower order model has provided quite accurate results. 

When a linear variation of the streamfunctions is assumed for a tri- 
angle as shown in figure 14, the streamfunction at any point (x, y) can 
be written as 


>Kx,y) = [B^ B^] 


B.tjj. 

— 1—1 


(A2) 


In equation (A2) the elements of row rector are first order poly- 

nomials in terms of x and y. The derivative of the streamfunction ip in 
X and y directions can be calculated as follows: 
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[Bi B- B_ ] 
*■ lx 2x Sx-* 


B. \p. 

— IX— 1 


(A3) 


and 


= [B B B ] 
3y '■ ly 2y 5y‘ 


= B. ip. 
— ly — 1 


(A4) 


The integral in equation (Al) can now be evaluated in a series form as 


4’(>i') =7 I B. + b!” B. ] ij;. a. 

/ , —1 I— IX —IX — ly — ly j —1 1 

i=l 


(A5) 


where the subscript i indicates summation of respective terms for each of 
the finite element triangles. 

After applying the principle of minimum potential energy, the solu- 
tion of Laplace's equation for the defined finite element grid can be ob- 
tained from the following system of linear algebraic equations: 

n 

[ 0 ] = 


Z 

i=l 


A. f B^ B. + B^ B. 

1 [ —IX —IX — ly — ly 




A \p 


(A6) 


The matrix A resulting from equation (A6) is a singular matrix. The 
boundary conditions around the airfoil and at the free stream boundaries 
have to be applied to obtain a unique solution to the problem. Once the 
stream functions at the nodes of the triangles are calculated, a set of 
streamlines passing above and under the oscillating airfoil can be drawn 
as the contour lines for the calculated values of stream functions (fig- 
ure 12) . 


Boundary Conditions 

The boundary conditions shown in equation (3) can be described in two 
steps . 
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Boundary conditions on the airfoil . - According to the potential flow 
theory, the velocity vectors on the airfoil boundary are tangent to the 
airfoil surface. In other words, the boundary of the airfoil forms a 
streamline. This boundary condition can be imposed on the discretized 
system by introducing the following relation: 

^ - jjj = 0 (A7) 

p q 

The constraining relation (A7) simply implies that stream functions at 
neighboring nodes are equal. 

Boundary conditions for free stream boundary . - The second type of 
boundary conditions are those defining free stream conditions around the 
considered flow field (outer boundaries). The flow around the airfoil 
can be represented for various angles of attack by specifying the stream 
functions at the nodes of the rectangular outer boundary, in desired 
positions. In a computational sense, this can be considered equivalent 
to rotating the entire rectangular frame in figure 2, instead of oscil- 
lating the airfoil. Thus, it is possible to use the same matrix A in 
equation (A6) throughout the analysis and save considerable computer 
time. 


Solution of the System of Equations 

The boundary conditions corresponding to the system of equations 
shown in equation (A6) can be added, as an additional constraint, in 
matrix form as follows: 


= 0 


(A8) 


on the airfoil and 


1= d (A9) 

on the free stream boundaries. The total system of equations for the po- 
tential flow analysis can now be written as follows: 



(AlO) 
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In potential flow analysis, the A, £ and ^ matrices are constant while 
the matrix ^ changes with the position of the airfoil. If the A and £ 
matrices are combined, equation (AlO) can be expressed as follows: 


* 

A 

g 


★ 

i 


h 

t 

g. 

0 




0 


where 


(All) 



c 





(A12) 


Here ^ is a column matrix consisting of Lagrangian multipliers corre- 
sponding to the boundary conditions on the airfoil surface. In this case 
* ^ 
matrix A corresponds to the free stream flow, is not singular, and can 
be inverted to obtain the stream functions for the free stream flow. 


= E ■ 3. aJ'^ g^E (A13) 

where 

*-l 

E = A h 

*-l 

and E “ ^ S. 

* 

Since matrices A and h are constant for the potential flow, at each time 
interval only the matrices g and £ have to be calculated. In cases where 
the entire triangular grid was rotated to represent the motion of the air- 
foil, the matrix ^ becomes a constant matrix, while the £ and d matrices 
are calculated at each time interval. 

The effects of different types of obstacles in the flow field can be 
determined from the solution routine by defining a different ^ matrix and 
*-l 

using the same A matrix. Here again the b matrix represents the geometry 
of the obstacle and the A matrix represents the free stream flow. 


(A14) 

(A15) 
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APPENDIX B 


FINITE ELEMENT ANALYSIS OF UNSTEADY 
FLOW AROUND AN OSCILLATING AIRFOIL 


Mathematical Description of the Problem 


Governing differential equations . - To analyze unsteady flow around 
an oscillating airfoil, one has to solve the Navier-Stokes equations for 
incompressible viscous fluid, which can be written as follows: 


3u 

3t 

3u 
“ 3x 

3u 
"" 3y 

= ix 

P 

1 ^ 
p 3x 

3v 

3t 

3v 

+ u 

3x 

3y 

= Iy 

p 

1 ^ 
" p 3y 




The required additional condition of continuity can be written as 


(Bl) 


(B2) 


3u 3v 
3x * 3y 


0 


(B3) 


The boundary conditions on the airfoil can be defined as u = v = 0. The 
velocity distribution on the outer rectangular boundary can be prescribed 
from the free stream conditions. 


Equations (Bl) and (B2) can now be expressed in a convenient form by 
introducing an additional variable to define vorticity as 


3v _ ^ 
3x 3y 


(B4) 


Equations (Bl) and (B2) can be combined as the Helmholtz equation 


3(i) 3 I 3(0 \ 3 ( 3(0 

+ -r— U(0 - V — + V(0 - V — 

3t 3x I 3x / 3y \ 3y 


= 0 


(B5) 


Equation (B5) is also known as the "vorticity transport" equation. 
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To satisfy the equation of continuity another variable will be 
introduced by defining the stream function as 


8y 3x 

The equation of continuity is then automatically satisfied. From equation 
(B4) the vorticity can be related to the stream function as follows: 

0) = - (B7) 

The values of the stream functions are known for the free stream boundary 
and are made equal to each other on the airfoil boundary. 

From the above formulation the finite element method can be applied for 
the solution of incompressible flow problems. Examples of such applications 
are these: 

a) potential flow; can be analyzed with equation (B7) for vanishing 
vorticity terms (Appendix A) 

b) steady state flow; can be analyzed from equations (B5) and (B7) 
by numerical integration of equation (B5) until a stable solution 
is obtained 

c} unsteady flow; requires the simultaneous solution of equations 
(B5) and (B7) and numerical integration of equation (B5). 

Once the velocity distribution has been determined, the corresponding pres- 
sure distribution may also be calculated from the following expression: 


2 2 

- 

2 2 

ax ay 


2 2 2 2 2 
ja V ^ 2 9 uv + 3 V 


ax 


ax3y 


ay 


(B8) 


Method of solution . - The mathematical analysis of unsteady flow prob- 
lems can be carried out in the following manner: 


a) determine initial vorticity distribution 

b) calculate stream functions from equation (B7) 

c) determine velocity distribution in x- and y- directions from 
equation (B6) 
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d) integrate equation (B5) with respect to time to obtain the 
vorticity distribution at time t + At 

e) repeat steps b, c, and d for each time increment. 

The initial vorticity distribution is usually unknown. Therefore start- 
ing from an arbitrary vorticity distribution for the initial stationary air- 
foil position, the procedure is repeated until a steady state solution is 
obtained. 

A numerical procedure based on the finite element method, was applied 
for each one of the above steps. 


Analysis of the Vorticity Transport Equation 

A finite element representation of an oscillating airfoil and the 
surrounding air flow is shown in figure 3 by a triangular mesh system. 
Velocity and vorticity distributions around the airfoil are defined at the 
nodes for each triangular element. The numerical method can be outlined in 
the following steps. 

The finite element . -A basic finite element model was chosen initially 
to examine the method of solution. It was assumed that vorticity and stream 
function values vary linearly over the area of each triangular element. 

Then, for a typical triangular element, shown in figure 15, the vorticity and 
the stream function at any point on the triangle are defined by 


and 


(u = Bj B3 


B.co. 

—1—1 


(B9) 



(BIO) 


The solution of Poisson's equation . - Knowing the vorticity distribu- 
tion around the airfoil, one can determine the stream functions from the fi- 
nite element solution of Poisson's equation. The solution of this equation 
can be obtained as described in Appendix A. The quadratic function to be 
minimized is obtained from the variational formulation of equation (B7) as 
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(Bll) 


(B12) 


The solution of equation (B12) provides the stream functions at the nodal 
points, and thus the velocity distribution around the airfoil. The next 
step is the numerical integration of equation (B5). 

Variational formulation of the vorticity transport equation . - The ex- 
act variational formulation of the vorticity transport equation does not ex- 
ist due to the nonlinearities in the equation; containing both velocities 
and derivatives of vorticities. Although the linear terms in equation (B5) 
can be written in variational form quite readily, one has to treat the velo- 
cities u and V as functions of the vorticity distribution to. A perturbation 
technique for treating such nonlinearities has been applied by Ecer [ref. 21] 
to finite element analysis of buckling problems in solid mechanics and has 
produced satisfactory results. For the variational formulation of the vor- 
ticity transport equation, a similar technique can be applied. By using a 
Taylor series expansion of the velocities in terms of the vorticity distribu- 
tion, at a particular time t^, the velocity distribution becomes 



The problem is thus reduced to calculating derivatives of the velocity dis- 
tribution in terms of vorticities. This can be done readily from the finite 
element formulation of Poisson's equation. ITie equations of continuity and 
vorticity can be written as 


V* u = 0 


from which 


3u 



V X u = 0 ) (B14] 


3u 

V X - = 1 (B15) 
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Using equation (B15) calculation of the first derivative of the velocity 
vector can be obtained from the solution of Poisson's equation for a unit 
vorticity, at each of the nodes of the triangular grid. Higher order 
derivatives can also be defined in a similar manner, again from the solu- 
tion of Poisson's equation. As described in Appendix A, this procedure 
does not require excessive amount of computations, since the solution of 
the equations for free stream conditions (A13) is already stored in the 
computer. Only the first term in the Taylor series was initially kept, 
since the velocity distribution is not very sensitive to the changes in 
vorticity distribution. The integral functional to be minimized in this 
case can be written as follows: 



For a triangular element, with a linear variation of vorticity and stream 
functions, equation (B16) can be discretized in finite element form as 
follows : 



The vorticity transport equation 
tion of the system of difference 


B^ B 

-y -y 


+ u B^ B + V B^ B WA CO 
— —X — ^ 


can then be solved by numerical integra- 
equations in (B17) . 


(B17) 


Boundary conditions . - The boundary conditions to be satisfied by the 
finite element model are these: 


^ - prescribed on the rectangular outer boundary as shown in figure 3 
ip - constant on the airfoil boundary 

(ji) = 0 on the rectangular outer boundary and on the boundary of the 
airfoil . 


A Lagrange multiplier technique was used in the solution of the un- 
steady flow equations to satisfy the above boundary conditions. 

A different grid was prepared for the finite element analysis as 
shown in figure 3. Elements of smaller size were used for the inner region 
of the rectangular grid, where the airfoil was pitching and plunging. The 
shape of the airfoil was defined by specifying coordinates of a series of 
points where the airfoil boundary conditions are satisfied. The boundary 
conditions for stream functions around the airfoil can be satisfied by 
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requiring that all points on the airfoil boundary have equal values. The 
boundary conditions used to specify vorticities on the airfoil and on the 
rectangular outer boundary can be satisfied by setting the vorticity values 
at these points equal to zero. All these boundary conditions are included 
in the solution of Poisson's equation and the vorticity transport equation 
as supplementary constraint equations. As the position of the airfoil 
changes, coordinates of new points are defined to describe the airfoil; 
thus only the constraint equations need to be changed at different time in- 
crements . 

The direction of the velocity vectors on the boundary of a moving air- 
foil is a function of the motion of the airfoil. This problem can thus be 
treated by using a higher order finite element model, having velocities as 
unknowns at each node. In such cases the direction of the velocity vec- 
tors can be defined at each node, considering also the motion of the air- 
foil. As the numerical integration interval becomes smaller, however, both 
solutions converge. 

Nimierical integration of the vorticity transport equation . - The numer- 
ical solution of the difference equation (B17) can be obtained by applying 
well-known numerical techniques [ref. 22] . Fromm [ref. 3] has experimented 
with various difference techniques and has studied the stability of these 
equations with respect to various parameters. Roache and Mueller [ref. 6] 
applied a similar method for the solution of the vorticity transport equation, 
as obtained from a finite difference approximation. However, none of the 
above attempts yielded the clear and effective representation of governing 
time-dependent equations as obtained from the finite element method. Time- 
dependent boundary conditions can be imposed on equation (B17) as an 
additional constraint, as explained above. ITie change of airfoil boundary 
conditions with respect to time can be introduced by specifying that the 
vorticities are equal to zero for the new position of the airfoil at time 
t = t^ + At. Using an Euler method of integration, the vorticity distribu- 
tion at time t + At can be written as 


(0 


t + At 



(B18) 


or writing equation (B18) in matrix form as , 

When the vorticity transport equation is solved from the finite element 
formulation to obtain 
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3b) 


at t = t 

o 

the problem reduces to numerical integration of a nonlinear discrete system 
with respect to time. The nonlinear system is defined for the nth time step 
by 


dm" 

dt 


A 0) 


(B20) 


and under the boundary conditions 


dt^ 1 


dt 


At 


n+1 



CB21) 


Standard numerical techniques were applied for solving the above system of 
difference equations (B20) and (B21). Formulations based on a finite dif- 
ference approximation become intractable as a higher order integration tech- 
nique is used for the solution of the time-dependent problem. In the pre- 
sented approach, based on a finite element formulation, time and space varia- 
bles are separated. ITie resulting difference equations can be integrated 
by any of the recognized numerical integration techniques, where the ele- 
ments of an error analysis are well defined. 


Calculation of pressures around the airfoil . - Finally, the pressures 
around the airfoil can be calculated from the finite element solution of 
equation (B8). To increase the accuracy of the solution, a higher order 
finite element model can be employed. An improved finite element model 
requires a higher order approximation of the stream functions, while allow- 
ing vorticities to vary linearly over the triangular element. The stream 
functions and their derivatives at the nodes can be defined as 

"^1 33T" 3^ ^2 JjT 3^ ’^3 3x 

In the above example the number of equations in matrix equation (B7) will 
be increased three times. However, the method of solution remains the 
same. The distribution of velocity is then assumed to be linear over the 
triangular element as follows; 


(B22) 
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Equation (B8) can be written in variational form as follows: 


B u . 


VI 


(B23) 


$ 



+ 




(B24) 


Assuming also a linear variation of pressure over the triangular element, 
the problem again reduces to the solution of Poisson's equation, where the 
right hand side of the resulting matrix equation can be expressed by a 
column vector as follows: 


Q - P 


u^. B^ B u, . 
•^1 —X —X 


2u^. B^ B u . + u^. B^ B u . 
—X ^ -^1 -^1 — y -^1 


(B25) 


After the pressure distributions around the airfoil are calculated, the 
lift and drag coefficients for the airfoil's unsteady pitching and plung- 
ing motion can be computed. VVhen computed results are tabulated, the un- 
steady aerodynamic characteristics of airfoils for various flight condi- 
tions are readily evaluated. 


Accuracy and Stability of Solutions of the Vorticity 
Transport Equation 

The numerical errors involved in the mathematical analysis of the vorticity 
transport equation can be summarized as computational errors and dis- 
cretization errors. Both kinds of errors are closely related to the 
stability of the governing differential equations and can be analyzed in 
two parts : 


. solution of the system of linear algebraic equations resulting 
from the finite element representation of Laplace's equation 

finite element formulation and numerical integration of the 
vorticity transport equation. 

Errors involved in the solution of Laplace's equation . - The 
Lagrange multiplier technique used in the solution of Laplace's equation 

aids in understanding the errors involved in the analysis. Matrix A* in 
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equation (AiO) represents Laplace's equation over a rectangular grid 
with specified boundary conditions. TTiis is equivalent to the finite 
element representation of the plane stress problem of a rectangular section. 
The accuracy of the solution has been discussed by various authors [ref. 23, 

2 

24]. The conditioning number for such a matrix is proportional to n , 
where n is the number of elements in the larger dimension of a rectangular 
grid. Consequently, for the grid used in the presented analysis, as shown 
in figure 3, satisfactory answers were obtained. 

Errors involved in the numerical integration of the vorticity transport 
equation . - The numerical integration of the vorticity transport equa- 
tion involved both computational and discretization errors. The accu- 
racy of the numerical integration can be checked by using different nu- 
merical integration techniques and different integration step sizes. 

The results in figure 16 indicate that results obtained from the numer- 
ical method are stable for the chosen particular time step sizes, 
depending on the position of a point with respect to the airfoil boundary. 

In order to understand this behavior, an error analysis can be made on a- 
typical triangular element of the grid in figure 3. 

The Euler method was chosen in illustrating the stability of the nu- 
merical results because of its simplicity. The instability conditions for 
a steady flow case are analyzed by neglecting the effects of generalized 
unsteady forces on the airfoil boundary, since they decrease rapidly as 
the solution converges to a steady state. The vorticity transport equa- 
tion can be written in matrix form as 

9w 

■^ = A w (B26) 

From which, using a forward difference relationship one can write 


“n + 1 = (1 A (B27) 

Following Richtmeyer [ref. 22] the stability of equation (B17) can be 
determined by checking whether the following expression has a finite value: 

lim (^ = (1 + A At)’^ ^ CB28) 

n — 
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Since matrix A converges to a constant matrix, the changes in this matrix can 
be neglected for the steady state problem, llie stability of equation (B17) 
can then be determined by checking whether the largest eigenvalue of 


or 


(^ + A At) 1 


(B29) 


At <_ _2 
6 

max 

^max estimated from the finite element model by obtaining an upper 

bound for the largest eigenvalue of matrix A. Since the system of equa- 
tions is assembled from triangular elements [ref. 24] 


^max— ^max (B30) 

where is the maximum value of the eigenvalues of each triangular 

element. The problem then reduces to estimating the largest eigenvalue 
of a typical triangle from the following equation 


1 max = 


0). A. 0). 
— 1 —1 —1 

t 

0). U). 

—1 —1 


(B31) 


For a triangular element, with a linear variation of vorticity and velocity 
distribution, the largest eigenvalue can be calculated by substituting 


0 ) 


0) 

-o 



(B32) 


into the vorticity transport equation. Equation [B17) can be written in 
matrix form for each triangular element as follows: 
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For a typical triangular element in figure 3, with base b and height h = 2b, 
the magnitude of the maximum eigenvalue of an element will depend on the 

kinematic viscosity v and the velocity distribution u and v. The stability 
of the solution for the vorticity distribution can be estimated from the 
largest eigenvalue of the system of equation CB33) for several typical cases . 


For a triangular element in the free stream, where u is constant and v = 0, 
the eigenvalues do not depend on the velocity distribution. The largest 
eigenvalue then becomes 
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CB34) 


As can be seen from equation (B34) , when the viscosity terms are neglected, 
the numerical integration of the inviscid flow equations becomes infinitely 
stable for the free stream flow. The importance of the velocity terms, 
however, is apparent at the triangles on the boundary of an obstacle. 
Assuming a triangle with a base on the airfoil boundary = u^ = 0) and 

neglecting the vertical velocities (v^ = V 2 = v^ = 0 ) , the largest eigen- 
value of equation (B33) can be written as 
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(B35) 








The stability of equation (B26) depends on the real part of the y. in equa- 
tion (B35). If the velocity term U 2 satisfies equation (B36), ^ 

u^ < 19.5 ^ (B36) 

then is a real number and the time step of the numerical integration can 

be estimated from equation (B36) as an upper bound. However, if equation 
(B36) is not satisfied, then the numerical integration will have an oscil- 
latory convergence. Also, equation (B36) shows that at the boundary layer 
region, where the velocity distribution varies parabolically starting from 
the airfoil boundary, stable results can be obtained by using smaller tri- 
angles and decreasing u^ for the triangles at the boundary of the airfoil. 
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Reynolds 

Number 

Angle of Attack 

Number of Integration 
Steps Required to 
Obtain a Steady State 

C.P.U. 

Time 

1000 

4 

140 

9.1 min 

1000 

8 

140 

10. 1 min 

105 

4 

75 

5 . 0 min 

105 

8 

76 

6.5 min 


Table 1. Computer Time Required for the Solution 
of Vorticity Transport Equation 












Fig. la. Flow chart of the computer program for the mathematical 
analysis and visualization of flow around oscillating airfoils 



Potential Flow 



Unsteady Viscous Flow 


Fig. lb. Flow chart of the program for the solution of the flow 
around an oscillating airfoil from potential flow and the 
vorticity transport equation 
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velocity 

distribution 


numerical integration 
of velocity distribution 
to determine streamlines, 
streaklines and pathlines 


motion of 
fluid particles 


computer graphics routines 
for generating illustrations 
and movies 


computer- drawn 


Fig. Ic. Flow chart of the program for generating computer-drawn 
movies from the output of the flow analysis 
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Fig. 2a. Finite Element Gridwork 





00 



Fig. 2b. Finite Element Gridwork for Area A 



Fig. 2c. Finite Element Gridwork for Area B 
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Fig.' 3. Complete Finite Element Gridwork for the Unsteady Viscous Flow Problem 



iig. 4a. Sample frames from a computer- drawn movie 
showing the flow of air particles around a pitching 
and plunging NACA 0012 airfoil. The mathematical 
model used for the preparation of the movie was 
based on potential flow theory. Each frame 
represents location of air particles at 1/24 
second time intervals. 







Fig. 4b. Sample computer-drawn illustration showing details 
of vorticity distribution around NACA 0012 airfoil 





Fig. 5a . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10'^, a = 0°) A vortex approaches and strikes the 
airfoil, starting from a position of 47% chordwidth ahead 
of the airfoil and separates into two vortices. The 
maximum vorticity on the vortex is initially equal to unity. 
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t = 3.2 sec. 


t = 4.0 sec. 


Fig. 5b . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a = 0°) . The vorticity transport equations were 
integrated with respect to time to obtain steady state 
conditions; starting from an arbitrary vorticity distribution 
i.e. two vortices, one below and one above the airfoil. 
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t = 4.8 sec. 


t = 5 . 2 sec . 




t = 5.6 sec. 


t = 6.4 sec. 




t = 7.2 sec. 


t = 8.0 sec. 


Fig. 5c . 2 Vorticity distribution around a NACA 0012 airfoil. 

(Re = 10 , a = 0°). Two vortices, advancing along the airfoil, 
converge to a steady state condition. Last frame shows the 
vortices at trailing and leading edges of the airfoil for 
steady state conditions. 
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t = 0.064 sec. t = 0.080 sec. 


Fig. 6a. g Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a =0°). A unit vortex, as in Fig. 5a, advances 
and strikes the airfoil. Compared to conditions at lower 
Reynolds numbers, the maximum vorticities are travelling 
closer to the airfoil. 


46 









t = 0.160 sec. 


t = 0.176 sec. 


Fig. 6b. g Vorticity distribution around a NACA 0012 airfoil. 

[Re = 10 , a = 0°). As two vortices advance along the airfoil, 
the magnitude of the vorticities at the leading edge of the 
airfoil increase. The vorticity distribution is relatively 
unstable at the leading edge. 
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t = 0.256 sec. 


t = 0.272 sec. 


Fig. 6c. g Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10^, a = 0°) . As the vorticity distribution around the 
airfoil reaches a steady state at the trailing edge , at the 
leading edge it is still unsteady. The last frame, for 
t = 0.272 sec., is drawn at a larger scale and shows that 
the vortex rings are distributed along the airfoil rather 
than concentrated at the trailing edge, as in Fig. 5c. 


48 










t = 0.1 sec. 


t = 0.8 sec. 




t = 1 . 6 sec . 


t = 2.4 sec. 




t = 3.2 sec. 


t = 4.0 sec. 


Fig. 7a . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a = 4°). A unit vortex, similar to the one in 
Fig. 5a, strikes the pitched airfoil. 
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t = 8.8 sec. 


Fig. 7b . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = lO"’, a = 4°). The initial unit vortex separates into 
two asymmetric vortices; upper vortex moves away from the 
airfoil while lower one reaches a maximum vorticity near the 
airfoil . 
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t = 12.8 sec. t = 13.6 sec. 


Fig. 7c . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a = 4°) . The numerical integration procedure was 
continued until a steady state distribution was obtained. 
Compared to Fig. 5c, these vortices are more evenly 
distributed along the airfoil and closer to the leading 
edge. Beneath the airfoil the vortices adhere closer 
to the surface and are larger in magnitude. 
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t = 0.040 sec. 


t = 0.048 sec. 


Fig. 8a. ^ Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a = 4°). A unit vortex, similar to the one in 
Fig. 6a, strikes the pitched airfoil. 
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t = 0.088 sec. 


t = 0.096 sec. 


Fig. 8b. ^ Vorticity distribution around a NACA 0012- airfoil . 
(Re = 10^, a = 4°). The vortex separates and advances along 
the airfoil. Compared to Fig. 6b, these two vortices are not 
symmetric. Compared to Fig. 7b, the upper vortex is closer 
to the airfoil and larger in magnitude. 


53 










t = 0.104 sec. 


t = 0.112 sec. 




t = 0.120 sec. 


t = 0.128 sec. 




t = 0.136 sec. 


t = 0.144 sec. 


Fig. 8c. g Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10^, a = 4°) . The vorticity distribution around the 
airfoil reaches a steady state. Compared to Fig. 7c, the 
vortices are much closer to the airfoil and exhibit an 
unsteady behavior at the leading edge. 
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t = 3.2 sec. 


t = 4.0 sec. 


Fig. 9a . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10 , a = 8°). A unit vortex, starting from the same 
position as in Fig. 5a and 7a, approaches the airfoil. 









t = 8.0 sec. 


t = 8.8 sec. 


Fig. 9b. Vorticity distribution around a NACA 0012 airfoil. 
(Re = lO"^, a = 8°). The vortex strikes the airfoil and 
separates into two vortices. Compared to conditions shown 
in Fig. 7b, the upper vortex remains closer to the airfoil. 
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Fig. 9c . 2 Vorticity distribution around a NACA 0012 airfoil. 
(Re = lO"^, a = 8°]. The vortices advance along the blade 
and reach a steady state at almost the same time as the case 
shown in Fig. 7c. The vortex generated in the unsteady 
region at the leading edge of the airfoil, as seen in Fig. 

5c and 7c, moves along the lower edge in this case. 
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t = 0.048 sec. 


t = 0.056 sec. 


Fig. 10a. Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10^, a = 8°) . A unit vortex, starting from the same 
position as in Fig. 6a and 8a, approaches the airfoil. 
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t = 0.064 sec. 


t = 0.072 sec. 






Fig. 10b. Vorticity distribution around a NACA 0012 airfoil. 
(Re = 10^, a = 8°) . The vorticity distribution at the leading 
edge of the airfoil is increased. 
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t = 0.144 sec. 


t = 0.152 sec. 


Fig. 10c. Vorticity distribution around a NACA 0012 airfoil. 

(Re = 10^, a = 8°). Steady state is obtained along the airfoil. 
Compared to Fig. 9c, the upper vortex is nearer to the trailing 
edge, while the lower vortex is lagging. The vorticity 
distribution is increased at the leading edge. 
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t = 0.012 sec. 


t = 0.024 sec. 




t = 0.036 sec. 


t = 0.048 sec. 




t = 0.060 sec. 


t = 0.072 sec. 


Fig. 11a. Vorticity distribution around a pitching NACA 0012 
airfoil. (Re = 10^). The airfoil is pitching upwards starting 
from a horizontal position. The initial vortex has the same 
magnitude as in the previous cases, however, the vorticity 
distribution was plotted using a larger scale. 
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t = 0.132 sec. 


t = 0.144 sec. 


Fig. 11b. Vorticity distribution around a pitching NACA 0012 
airfoil. (Re = 10^). As the vortex strikes the airfoil a 
series of vortices are generated. These vortices travel 
along the blade in a similar fashion to the ones in Fig. 10b. 
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t = 0.204 sec. 


t = 0.216 sec. 


Fig. 11c. Vorticity distribution around a pitching NACA 0012 
airfoil. (Re = 10^). The numerical integration of the vorticity 
transport equations was advanced in time for the case shown in 
Fig. 11a and 11b. The vorticity distribution was plotted using 
a larger scale. As the blade continues pitching, new vortices 
are generated at the leading edge of the blade. 
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Fig. 12. Streamlines around a NACA 0012 airfoil. The 
streamlines were plotted from values of streamfunctions at 
each step of the numerical integration of the vorticity 
transport equations. For the sample cases shown in 
Fig. 5 - 1, streamlines are illustrated above at several 

time increments. 
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Fig. 13. Streaklines around a NACA 0012 airfoil. The velocity 
distribution around the airfoil was obtained from the solution 
of vorticity transport equations, as a function of time. The 
positions of a series of air particles were calculated by 
numerical integration, the position of several particles 
which started from the same points were plotted. 
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Vorticity, sec. Vorticity, sec. 



Time, sec. 

Fig. 16a. Accuracy and stability of numerical integration of the 
vorticity transport equation for a point (A) upstream of the airfoil. 

(a = 0, Re = 10^). 
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Fig. 16b. Accuracy and stability of numerical integration of the 
vorticity transport equation for a point (B) in the vicinity of the 

airfoil. (a = 0, Re = 10^). 
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